subroutine set_vbox(lunit,t_array0,myrank)
  use param
  use box_param,only:bwidth
  use box_param,only:bound
  use box_param,only:cell_width
  use box_param,only:vbox_map
  use box_param,only:vbox,n_vbox
  use react_data_param,only:react_type
  use react_data_param,only:n_react_type
  implicit none

  interface
     function reallocReal(p,n)
       use param
       real(dbl),pointer,dimension(:) :: reallocReal,p
       integer,intent(in) :: n
     end function reallocReal

     function reallocInt(p,n)
       integer,pointer,dimension(:) :: reallocInt,p
       integer,intent(in) :: n
     end function reallocInt
  end interface

  integer,intent(in) :: lunit,t_array0(8),myrank

  integer :: i,j,k,n,p
  integer :: ni,nj,nk,tx,ty,tz
  integer :: ierr
  integer :: nx,ny,nz,nn
  real(dbl) :: maxr
  character(len=8) :: d_time_str

  ! If the number of the react_type is zero,
  ! it doesn't need to set the virtual box.
  ! If the smallest size of the virtual box
  ! is smaller than maximum radius of interaction,
  ! the size of the cell is replaced by the maxr.

  ! Determine nx,ny,nz for different cases
  if (n_react_type.eq.0) then
     ! No need to build up the neighbor list
     nx=1
     ny=1
     nz=1
     ! cell_width is set as the width of box
     maxr=bwidth(1)
     if (bwidth(2).gt.bwidth(1)) maxr=bwidth(2)
     if (bwidth(3).gt.bwidth(2)) maxr=bwidth(3)
     cell_width=maxr
  else
     ! Find the maximum capture radius.
     maxr=react_type(1)%react_rd(1)
     do i=2,react_type(1)%n_data
        if (maxr.lt.react_type(1)%react_rd(i)) &
             maxr=react_type(1)%react_rd(i)
     end do
     do i=2,n_react_type
        do j=1,react_type(i)%n_data
           if (maxr.lt.react_type(i)%react_rd(j)) &
                maxr=react_type(i)%react_rd(j)
        end do
     end do
     ! maxr is set a little bit larger than maximum
     ! capture radius.
     maxr=maxr+0.1*maxr
     ! Determine the size of the virtual cell
     if (cell_width.lt.maxr) cell_width=maxr
     nx=ceiling(bwidth(1)/cell_width)
     ny=ceiling(bwidth(2)/cell_width)
     nz=ceiling(bwidth(3)/cell_width)
  end if

  ! Mapping function
  ! (ix,iy,iz) <-> 1d index
  ! At the boundary cell, vbox_map(i,j,k) is equal to 0.
  allocate(vbox_map(0:nx+1,0:ny+1,0:nz+1),stat=ierr)
  if (ierr.ne.0) call abnormal_stop("set_vbox",&
       "Memory allocation to vbox_map has failed.")
  do k=0,nz+1
     do j=0,ny+1
        do i=0,nx+1
           vbox_map(i,j,k)=0
        end do
     end do
  end do

  ! Initialize vbox
  allocate(vbox(nx*ny*nz),stat=ierr)
  if (ierr.ne.0) call abnormal_stop("set_vbox",&
       "Memory allocation to vbox has failed.")

  n_vbox=0

  do k=1,nz
     do j=1,ny
        do i=1,nx
           n_vbox=n_vbox+1
           ! save the mapping index
           vbox_map(i,j,k)=n_vbox
           vbox(n_vbox)%pos(3)=k
           vbox(n_vbox)%pos(2)=j
           vbox(n_vbox)%pos(1)=i
           ! number of the objects belonging to this vbox
           vbox(n_vbox)%n_obj=0
           ! Indicate the index of the first object
           vbox(n_vbox)%first=0
           ! Save the neighbor objects to this box
           vbox(n_vbox)%n_neighbor=0
           nullify(vbox(n_vbox)%neighbor_list)
           nullify(vbox(n_vbox)%dx)
           nullify(vbox(n_vbox)%dy)
           nullify(vbox(n_vbox)%dz)
        end do
     end do
  end do

  ! Set the neighbor and treat the boundaries
  do n=1,n_vbox
     do k=vbox(n)%pos(3)-1,vbox(n)%pos(3)+1
        do j=vbox(n)%pos(2)-1,vbox(n)%pos(2)+1
           do i=vbox(n)%pos(1)-1,vbox(n)%pos(1)+1
              ! current index
              p=vbox_map(i,j,k)
              ! Index considering the boundary condition
              nk=k
              nj=j
              ni=i
              ! necessary when evaluate the distance 
              ! between objects
              tx=0
              ty=0
              tz=0
              ! If the neighbor cell is boundary,
              ! apply the boundary condition. (0/1)
              if (p.eq.0) then
                 ! Along the x-axis
                 if (i.eq.0) then
                    ! left boundary
                    if (bound(1).eq.1) ni=nx
                    if (bound(1).eq.1) tx=1
                 else if (i.eq.(nx+1)) then
                    ! right boundary
                    if (bound(1).eq.1) ni=1
                    if (bound(1).eq.1) tx=2
                 else
                    ! Inner box
                 end if

                 ! Along the y-axis
                 if (j.eq.0) then
                    ! left boundary
                    if (bound(2).eq.1) nj=ny
                    if (bound(2).eq.1) ty=1
                 else if (j.eq.(ny+1)) then
                    ! right boundary
                    if (bound(2).eq.1) nj=1
                    if (bound(2).eq.1) ty=2
                 else
                    ! Inner box
                 end if

                 ! Along the z-axis
                 if (k.eq.0) then
                    ! left boundary
                    if (bound(3).eq.1) nk=nz
                    if (bound(3).eq.1) tz=1
                 else if (k.eq.(nz+1)) then
                    ! right boundary
                    if (bound(3).eq.1) nk=1
                    if (bound(3).eq.1) tz=2
                 else
                    ! Inner cell
                 end if
              end if

              ! Skip for the free boundary 
              p=vbox_map(ni,nj,nk)
              if (p.eq.0) cycle

              vbox(n)%n_neighbor=vbox(n)%n_neighbor+1
              nn=vbox(n)%n_neighbor
              vbox(n)%neighbor_list=>&
                   reallocInt(vbox(n)%neighbor_list,nn)
              vbox(n)%neighbor_list(nn)=p
              vbox(n)%dx=>reallocReal(vbox(n)%dx,nn)
              vbox(n)%dy=>reallocReal(vbox(n)%dy,nn)
              vbox(n)%dz=>reallocReal(vbox(n)%dz,nn)

              ! Determine dx
              select case(tx)
              case(0)
                 vbox(n)%dx(nn)=0
              case(1)
                 vbox(n)%dx(nn)=-bwidth(1)
              case(2)
                 vbox(n)%dx(nn)=bwidth(1)
              end select

              ! Determine dy
              select case(ty)
              case(0)
                 vbox(n)%dy(nn)=0
              case(1)
                 vbox(n)%dy(nn)=-bwidth(2)
              case(2)
                 vbox(n)%dy(nn)=bwidth(2)
              end select

              ! Determine dz
              select case(tz)
              case(0)
                 vbox(n)%dz(nn)=0
              case(1)
                 vbox(n)%dz(nn)=-bwidth(3)
              case(2)
                 vbox(n)%dz(nn)=bwidth(3)
              end select
           end do
        end do
     end do
  end do

  if (n_react_type.eq.0) then
     if (n_vbox.ne.1) call abnormal_stop("set_vbox",&
          "n_vbox must be one if n_react_type=1.")
     vbox(1)%n_neighbor=0
  end if

  call print_d_time(t_array0,d_time_str)  
  if (myrank.eq.0) &
       write(*,'(a8,":SET_VBOX DONE.")') d_time_str(1:8)
  write(lunit,'(a8,":SET_VBOX DONE.")') d_time_str(1:8)

  return

end subroutine set_vbox
